%%% Matteo Ciccarelli and Benoit Mojon (c). Replicate results in Table 1-2

clear all;

%% load a panel of GDP for G7 countries
%data = xlsread('data section2 and section 4.xls');
data = xlsread('cpi_DATA.xls');
% ordering of the series
%series = {CPIQU2	CPIQUS	CPIQCA	CPIQUK	CPIQJP	CPIQDE	CPIQFR	CPIQIT	CPIQAT	CPIQBE	CPIQDK	CPIQFI	CPIQGR	CPIQIE	CPIQLU	CPIQPT	CPIQES	CPIQSE	CPIQNL	CPIQAU	CPIQNZ	CPIQNO	CPIQCH	CPIOECD

nco = 24;
Y = data(:,1:nco); 
n = size(Y,2);

%% compute annualized quarterly growth rates
dly = [log(Y(5:end,:))-log(Y(1:end-4,:))]*100;

%% use dly or dly1
x = dly(:,2:n-1);
[T,N] = size(x);

%%% principal components
OPTS.disp=0;
Mx = mean(x);
Wx = (std(x));

%%% center data
x_std = (x-kron(ones(T,1),Mx))*diag(1./Wx);

[v,d] = eigs(cov(x_std),1,'lm',OPTS);
F_pc = x_std*v*d^-.5;
chi_pc = x_std*v*v'*diag(Wx) + kron(ones(size(x,1),1),Mx);

%%% factor analysis
r = 2; %% number of factors

%[chi_zz,F_zz,C_zz,R_zz] = fa_zz(x,r); %% static em etimation
[chi_em,F_em,C_em,R_em] = FA_em(x,r); %% static zig-zag estimation
%[chi_dyn,F_dyn,A_dyn, C_dyn, Q_dyn, R_dyn] = dyn_em(x,r,0,r,1); %% dynamic factor analysis q = r = 1; p =1; s = 0
%[chi_dyn4,F_dyn4,A_dyn4, C_dyn4, Q_dyn4, R_dyn4] = dyn_em(x,r,1,r,1); %% dynamic factor analysis q = r = 1; p =2; s = 3

%% cross-sectional averages
avg_x = mean(x_std,2); avg_x = avg_x./std(avg_x);
chi_mean = avg_x*inv(avg_x'*avg_x)*avg_x'*x_std*diag(Wx) + kron(ones(size(x,1),1),Mx);
C_ave = (inv(avg_x'*avg_x)*avg_x'*x_std*diag(Wx))';
%%% compute sign to rotate estimated factors
%temp = corrcoef(avg_x,F_pc(:,1)); s_pc=sign(temp(2,1));
%temp = corrcoef(avg_x,F_zz(:,1)); s_zz=sign(temp(2,1));
%temp = corrcoef(avg_x,F_em(:,1)); s_em=sign(temp(2,1));
%temp = corrcoef(avg_x,F_dyn(:,1)); s_dyn=sign(temp(2,1));

% temp = corrcoef([avg_x,F_dyn]); s_dyn = sign(temp(1,2:end));  
% temp = corrcoef([avg_x,F_zz]); s_zz = sign(temp(1,2:end));  
temp = corrcoef([avg_x,F_em]); s_em = sign(temp(1,2:end));  
% temp = corrcoef(avg_x,F_pc(:,1)); s_pc=sign(temp(2,1));
if s_em(1,1)<0
    F_em(:,1) = s_em(1,1)*F_em(:,1);
end

    
    vari_em1=(C_em(:,1).^2)*cov(F_em(:,1))./diag(cov(x));
    vari_em2=(C_em(:,2).^2)*cov(F_em(:,2))./diag(cov(x));
%    disp('Variance decomposition')
%    disp([ vari_em1 vari_em2 ])
%    disp([ vari_em1 ])
 
% if r==1
%     F_dyn = s_dyn*F_dyn ;
%     F_em = s_em*F_em ;
%     vari_ave=(C_ave(:,1).^2)*cov(avg_x(:,1))./diag(cov(x));
%     vari_em=(C_em(:,1).^2)*cov(F_em(:,1))./diag(cov(x));
%     vari_dyn=(C_dyn(:,1).^2)*cov(F_dyn(:,1))./diag(cov(x));
%     disp([vari_ave vari_em vari_dyn])
% end 
%     
% if r==2
%     vari_dyn1=(C_dyn(:,1).^2)*cov(F_dyn(:,1))./diag(cov(x));
%     vari_dyn2=(C_dyn(:,2).^2)*cov(F_dyn(:,2))./diag(cov(x));
%     disp('Variance decomposition')
%     disp([ vari_dyn1 vari_dyn2 ])
% end
%%
x = dly(:,1:n-1);
[T,N] = size(x);
oecd_fact = (dly(:,nco)-kron(ones(T,1),mean(dly(:,nco))))*diag(1./std(dly(:,nco)));

regr = [ones(T,1) oecd_fact];

for i= 1:N
    par=inv(regr'*regr)*regr'*x(:,i);           %%% ols coeff
    u=x(:,i)-regr*par;                    %%% residuals
    m=size(regr); Ti=m(1); k=m(2);
    esu=u'*u/(Ti-k);
    r2_oecd(i,1)=1-(u'*u)/(center(x(:,i))'*center(x(:,i)));              %%% R2
    r2_oecd(i,2)=1-(1-r2_oecd(i,1))*(Ti-1)/(Ti-k);        %%% adjusted R2
end

regr_av = [ones(T,1) avg_x];

for i= 1:N
    par=inv(regr_av'*regr_av)*regr_av'*x(:,i);           %%% ols coeff
    u=x(:,i)-regr_av*par;                    %%% residuals
    m=size(regr_av); Ti=m(1); k=m(2);
    esu=u'*u/(Ti-k);
    r2_av(i,1)=1-(u'*u)/(center(x(:,i))'*center(x(:,i)));              %%% R2
    r2_av(i,2)=1-(1-r2_av(i,1))*(Ti-1)/(Ti-k);        %%% adjusted R2
end

regr_fa = [ones(T,1) F_em(:,1)];
par=inv(regr_fa'*regr_fa)*regr_fa'*x(:,1);           %%% ols coeff
u=x(:,1)-regr_fa*par;                    %%% residuals
m=size(regr_fa); Ti=m(1); k=m(2);
esu=u'*u/(Ti-k);
r2_eu1(1,1)=1-(u'*u)/(center(x(:,1))'*center(x(:,1)));              %%% R2
r2_eu1(1,2)=1-(1-r2_eu1(1,1))*(Ti-1)/(Ti-k);        %%% adjusted R2

regr_fa = [ones(T,1) F_em(:,2)];
par=inv(regr_fa'*regr_fa)*regr_fa'*x(:,1);           %%% ols coeff
u=x(:,1)-regr_fa*par;                    %%% residuals
m=size(regr_fa); Ti=m(1); k=m(2);
esu=u'*u/(Ti-k);
r2_eu2(1,1)=1-(u'*u)/(center(x(:,1))'*center(x(:,1)));              %%% R2
r2_eu2(1,2)=1-(1-r2_eu2(1,1))*(Ti-1)/(Ti-k);        %%% adjusted R2
% Table 1
sh_var_avoecd = [r2_av(:,1) r2_oecd(:,1)];
sh_var_fact = [vari_em1  vari_em2 ];

r2_em1 = zeros(size(sh_var_avoecd,1),1); 
r2_em1(1,1) = r2_eu1(1,1);
r2_em1(2:end,1) = vari_em1;

r2_em2 = zeros(size(sh_var_avoecd,1),1); 
r2_em2(1,1) = r2_eu2(1,1);
r2_em2(2:end,1) = vari_em2;

%%% this replicate table 1. it must then be sorted in ascending order along
%%% the first column
disp('table1')
sh_var_table1 = [r2_av(:,1) r2_oecd(:,1) r2_em1(:,1) r2_em2(:,1)]


%% table 2
% full sample
swit = 1; % 1 = G6; 2 = G8

if swit==1
    yy = [x(:,20) x(:,3) x(:,6) x(:,4) x(:,5) x(:,2)];
else
    yy = [x(:,20) x(:,3) x(:,6) x(:,4) x(:,5) x(:,2) x(:,7) x(:,8)];
end

Myy = mean(yy);
Wyy = (std(yy));

%%% center data
yy_std = (yy-kron(ones(T,1),Myy))*diag(1./Wyy);

avg_yy = mean(yy_std,2); avg_yy = avg_yy./std(avg_yy);

regr_av1 = [ones(T,1) avg_yy];

for i=1:N
    par=inv(regr_av1'*regr_av1)*regr_av1'*x(:,i);           %%% ols coeff
    u=x(:,i)-regr_av1*par;                    %%% residuals
    m=size(regr_av1); Ti=m(1); k=m(2);
    esu=u'*u/(Ti-k);
    r2_av1(i,1)=1-(u'*u)/(center(x(:,i))'*center(x(:,i)));              %%% R2
    r2_av1(i,2)=1-(1-r2_av1(i,1))*(Ti-1)/(Ti-k);        %%% adjusted R2
end

%% table 2
% 1975-2008
if swit==1
    yy = [x(57:T,20) x(57:T,3) x(57:T,6) x(57:T,4) x(57:T,5) x(57:T,2)];
else
    yy = [x(57:T,20) x(57:T,3) x(57:T,6) x(57:T,4) x(57:T,5) x(57:T,2) x(57:T,7) x(57:T,8)];
end

TT = size(yy,1);
Myy = mean(yy);
Wyy = (std(yy));

%%% center data
yy_std = (yy-kron(ones(TT,1),Myy))*diag(1./Wyy);

avg_yy = mean(yy_std,2); avg_yy = avg_yy./std(avg_yy);

regr_av1 = [ones(TT,1) avg_yy];

for i=1:N
    par=inv(regr_av1'*regr_av1)*regr_av1'*x(57:T,i);           %%% ols coeff
    u=x(57:T,i)-regr_av1*par;                    %%% residuals
    m=size(regr_av1); Ti=m(1); k=m(2);
    esu=u'*u/(Ti-k);
    r2_av2(i,1)=1-(u'*u)/(center(x(57:T,i))'*center(x(57:T,i)));              %%% R2
    r2_av2(i,2)=1-(1-r2_av1(i,1))*(Ti-1)/(Ti-k);        %%% adjusted R2
end

%% table 2
% 1983-2008
if swit==1
    yy = [x(89:T,20) x(89:T,3) x(89:T,6) x(89:T,4) x(89:T,5) x(89:T,2)];
else
    yy = [x(89:T,20) x(89:T,3) x(89:T,6) x(89:T,4) x(89:T,5) x(89:T,2) x(89:T,7) x(89:T,8)];
end

TT = size(yy,1);

Myy = mean(yy);
Wyy = (std(yy));

%%% center data
yy_std = (yy-kron(ones(TT,1),Myy))*diag(1./Wyy);

avg_yy = mean(yy_std,2); avg_yy = avg_yy./std(avg_yy);

regr_av1 = [ones(TT,1) avg_yy];

for i=1:N
    par=inv(regr_av1'*regr_av1)*regr_av1'*x(89:T,i);           %%% ols coeff
    u=x(89:T,i)-regr_av1*par;                    %%% residuals
    m=size(regr_av1); Ti=m(1); k=m(2);
    esu=u'*u/(Ti-k);
    r2_av3(i,1)=1-(u'*u)/(center(x(89:T,i))'*center(x(89:T,i)));              %%% R2
    r2_av3(i,2)=1-(1-r2_av1(i,1))*(Ti-1)/(Ti-k);        %%% adjusted R2
end

% show Table 2 for all countries with global inflAtion based only on g6 or
% g8. to get Table 2 select only the approriate countries
disp('table2')
tab2 = [r2_av1(:,1) r2_av2(:,1) r2_av3(:,1)]


%%%

% x = x-chi_mean;
% 
% [T,N] = size(x);
% 
% %%% principal components
% OPTS.disp=0;
% Mx = mean(x);
% Wx = (std(x));
% 
% %%% center data
% x_std = (x-kron(ones(T,1),Mx))*diag(1./Wx);
% 
% [v,d] = eigs(cov(x_std),1,'lm',OPTS);
% F_pc = x_std*v*d^-.5;
% chi_pc = x_std*v*v'*diag(Wx) + kron(ones(size(x,1),1),Mx);
% 
% %%% factor analysis
% r = 1; %% number of factors
% 
% %[chi_zz,F_zz,C_zz,R_zz] = fa_zz(x,r); %% static em etimation
% [chi_em,F_em,C_em,R_em] = FA_em(x,r); %% static zig-zag estimation
% %[chi_dyn,F_dyn,A_dyn, C_dyn, Q_dyn, R_dyn] = dyn_em(x,r,0,r,1); %% dynamic factor analysis q = r = 1; p =1; s = 0
% %[chi_dyn4,F_dyn4,A_dyn4, C_dyn4, Q_dyn4, R_dyn4] = dyn_em(x,r,1,r,1); %% dynamic factor analysis q = r = 1; p =2; s = 3
% 
% %% cross-sectional averages
% avg_x = mean(x_std,2); avg_x = avg_x./std(avg_x);
% chi_mean = avg_x*inv(avg_x'*avg_x)*avg_x'*x_std*diag(Wx) + kron(ones(size(x,1),1),Mx);
% C_ave = (inv(avg_x'*avg_x)*avg_x'*x_std*diag(Wx))';
% %%% compute sign to rotate estimated factors
% %temp = corrcoef(avg_x,F_pc(:,1)); s_pc=sign(temp(2,1));
% %temp = corrcoef(avg_x,F_zz(:,1)); s_zz=sign(temp(2,1));
% %temp = corrcoef(avg_x,F_em(:,1)); s_em=sign(temp(2,1));
% %temp = corrcoef(avg_x,F_dyn(:,1)); s_dyn=sign(temp(2,1));
% 
% % temp = corrcoef([avg_x,F_dyn]); s_dyn = sign(temp(1,2:end));  
% % temp = corrcoef([avg_x,F_zz]); s_zz = sign(temp(1,2:end));  
% temp = corrcoef([avg_x,F_em]); s_em = sign(temp(1,2:end));  
% % temp = corrcoef(avg_x,F_pc(:,1)); s_pc=sign(temp(2,1));
% if s_em(1,1)<0
%     F_em(:,1) = s_em(1,1)*F_em(:,1);
% end
% 
%     
%     vari_em1=(C_em(:,1).^2)*cov(F_em(:,1))./diag(cov(x));
%     disp('Variance decomposition')
%     disp([ vari_em1  ])
